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Abstract. We perform a multimode treatment of spin squeezing in- 
duced by interactions in atomic condensates, and we show that, at finite 
temperature, the maximum spin squeezing has a finite limit when the 
atom number A'^ — ^ oo at fixed density and interaction strength. To cal- 
culate the limit of the squeezing parameter for a spatially homogeneous 
■ system we perform a double expansion with two small parameters: 1/A'' 

^ I in the thermodynamic limit and the non-condensed fraction (Nnc) /N 

in the Bogoliubov limit. To test our analytical results beyond the Bo- 
goliubov approximation, and to perform numerical experiments, we use 
improved classical field simulations with a carefully chosen cut-off, such 
that the classical field model gives for the ideal Bose gas the correct 
non-condensed fraction in the Bose-condensed regime. 



1 Introduction 

A two-level atom can be described as an effective spin 1/2. Here, to describe an 
ensemble of atoms in two different internal states a and b, that are typically two 
hyperfine states, we use the picture of a "collective spin". This spin, of length N/2, is 
simply the sum of the effective spins 1/2 that describe the internal degrees of freedom 
of each atom. In the second quantized formalism the three hermitian spin components 
Sx, Sy and Sz are defined by: 

S+ s & + iS, = j d'r ild )itlt ) (1) 

S. - ^ (2) 

where the bosonic field operators 'ipa,b obey the usual commutation relations, Na = 
J (Pr il!\{r)tpa{T^) is the atom number in component a and the same for b. The spin 
operators are dimensionless and obey the commutation relations [Sx,Sy] = iSz and 
cyclic permutations. Physically Sz is the population difference between a and b states, 
while Sx and 5*1, describe one-body coherence between them. 
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Spin squeezing [T] is about creating quantum correlations, in such an ensemble of 
atoms, that can be useful for metrology. In particular spin squeezed states can be used 
to improve the accuracy of atomic clocks beyond the so called "standard quantum 
limit" that has been already reached in the most precise clocks |5] . The resulting gain 
for metrology is quantified by a spin squeezing parameter 



^ |(S)|2 ' ^""^ 

where N is the total atom number and AS'^ ^^^^ is the minimal variance of the collec- 
tive spin orthogonally to the direction of its mean value (S ) . The state is squeezed if 
and only if < 1 . As explained in |2 , in an atomic clock experiment using Ramsey 
population spectroscopy, ^ directly gives the reduction in the statistical fluctuations 
of the measured frequency cjab with respect to using uncorrelated atoms (for the same 
atom number N and the same Ramsey time T): 

^<i = ^^<r-^ (4) 



NT 

The parameter ^ in Eq.® is in fact the properly normalized ratio between the "noise" 
^S'.L^inin and the "signal" |(S)|. In experiments ^S'±,min is directly measured by 
measuring Sz after an appropriate state rotation and |(S)| is separately deduced from 
the Ramsey fringes contrast. 

Very recently experimental breakthroughs in spin squeezing have been achieved 
using either the interaction between atoms and light in an optical cavity [5] or atomic 
interactions in bimodal Bose-Einstein condensates [5], [7|- The ultimate limits of 
the different paths to spin squeezing are still objects of active studies [8, 9 10. ll| 12j . 
We address here the issue of non-zero temperature and of the influence of the non 
condensed fraction for spin squeezing schemes using Bose-Einstein condensates. 

We face the following physical problem: An interacting Bose gas, prepared at fi- 
nite temperature in the internal state a, is subjected to a sudden 7r/2 mixing pulse 
that puts each atom in a coherent superposition of two different internal states a and 
b. From this out of equilibrium state, with factorized spin and motional variables, 
quantum correlations and spin squeezing are created dynamically by the atomic in- 
teractions [I], Let us first sketch how this happens in a simple two- mode picture, 
i.e. assuming that all the atoms in a or 6 share the same wave function for their 
motional degrees of freedom. After the mixing pulse, the two condensates in a and 
b have a well defined relative phase, with a relative phase distribution whose width 
scales as and fluctuations in the relative particle number difference scaling 

as \/N. In a two-mode picture, the initial state can be expanded over Fock states 
\Na,N — Na) with Na particles in state a and N — Na particles in state b. Due to 
atom-atom interactions, each Fock state acquires a phase in the evolution that is 
proportional to Na — Ni, |lll3ll4j . This situation is completely equivalent to the evo- 
lution of a coherent state in a Kerr medium in optics. During the evolution, due to 
the different phase shifts of the different Fock states, the relative phase distribution 
starts to spread. At the same time, quantum fluctuations orthogonal to the mean 
spin direction get distorted and, before the relative phase distribution has sensibly 
spread, spin squeezing is created in the sample. Our aim is to include the two-mode 
quantum dynamics that we just described, and the effect of the thermally excited 
non-condensed modes within the same formalism. The thermal modes also provide 
a condensate phase spreading [H] , [H] , [HI , [IHl and are expected to affect the spin 
squeezing generated in the system at non-zero temperature [5] . For a review of spin 
squeezing and decoherence see also [TT] . 
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A central issue is the scaling of the squeezing as the system gets large, i.e. in 
the thermodynamic limit. Most studies are based on a two-mode description [1]. In 
this frame the squeezing parameter minimized over time tends to zero (infinite 
metrology gain) for A'^ — > oo as ^iin~iV-2/3. Although some studies beyond the 
two-mode theory were performed [4119120] they could not prove or disprove the two- 
mode scaling of spin squeezing in real condensates. Here we can go further. We find 
that for realistic atom numbers, the two-mode scaling ^^^^ ^ ]Sf-^/^ is meaningless at 
finite temperature and that the spin squeezing parameter ^^j^ at the thermodynamic 
limit has a finite non-zero value that we calculate explicitly. In this paper we present 
a detailed derivation of the results given in [5] and we present new improved classical 
field simulations, with a carefully chosen cut-off such that the classical field model 
gives for the ideal Bose gas the correct non-condensed fraction in the Bose-condensed 
regime. We also present results for the squeezing that would be measured by detecting 
only the condensed particles, which we call the "condensate squeezing" , and we show 
that it is much worse than the squeezing of the total field for reasons that we explain 
in the paper. 

In section [2] we formalize the problem and expose our approach to solve it. In 
section |3] we proceed with two numerical experiments. These experiments show (i) 
the existence of a non-zero thermodynamic limit for the squeezing parameter in con- 
trast with the predictions of the two-mode theory, and (ii) the universal scaling with 
the temperature of the squeezing in the thermodynamic and weakly interacting limit. 
Analytical calculations are performed in sectional By performing a double expansion 
of in terms of two small parameters, the inverse atom number 1/N controlling the 
thermodynamic limit and the non-condensed fraction controlling the weakly interact- 
ing limit, we obtain explicitly the minimal squeezing parameter that it is possible to 
achieve by this method as a function of the initial temperature and the interaction 
strength. A physical interpretation of the results is given in section [S] In that section 
we also show that the squeezing defined for the total field and the squeezing defined 
for the condensate mode only are very different and we give a physical explanation. 
We conclude in sectional 



2 The problem 

2.1 The Quantum Model 

We consider a spatially homogeneous system of A^ bosons in two internal states that 
interact with short range binary interactions. We take for simplicity identical in- 
teractions in components a and b and no crossed a-b interactions Q. The system is 
discretized on a cubic lattice of lattice spacing Z, with periodic boundary condition 
of period L along each direction x,y,z. For numerical convenience, L/l = rimax is 
an even integer. There are in total Af = V/dV lattice points, where V = is the 
system volume and dV = the unit cell volume. The Hamiltonian for one separate 
spin component, e.g. component a, reads 

^« = E + ^^t(r)v;t (r)V-,(r)^,(r) (5) 

k r 

In the kinetic energy term we have expanded the field operator over plane waves 



v-^ e 

V'a(r) = 2^ak-^ (6) 



^ In *^Rb atoms this may be done by spatial separation of the spin states [7| or by Feshbach 
tuning of the a-b scattering length [6]. 
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and Ok annihilates a particle of wave vector k belonging to the first Brillouin zone 
(FBZ) [— 7r/Z,7r/Zp of the lattice, so that along each direction ly, e ^{—L/{21), . . ., 
L/{21) — 1}. Since we consider a lattice model, the field operator here obeys the 
discrete bosonic commutation relations (j76p . The second term in (O represents atomic 
interactions modeled by a purely on-site interaction with a bare coupling constant on 
the lattice go. In practice, to recover the continuous space physics, I is taken to be 
smaller than both the healing length ^hcai and the thermal de Broglie wavelength AdB ■ 
In the weakly interacting regime, \a\ <^ Cheai, AdB, one can further take I ^ \a\ so that 
in the following we will identify go with the effective coupling constant g = Airh'^a/'m 
where a is the s-wave scattering length 0. 

Initially at t < 0, all the N atoms are in the internal state a in thermal equilibrium 
described by the canonical density operator 

P = ^e-^"^ (7) 

with /3 ~ l/{kBT) and T Tc where Tc is the transition temperature for Bose- 
Einstein condensation. 

At t an electromagnetic pulse mixes the states a and b. The pulse Hamiltonian 
acting during a time interval ^puisc is 

yp = ^T.'^yi^^l^b~i^Ua) (8) 

r 

In practice the timcscalc tpuiso is shorter than all the relevant timescales in the original 
Hamiltonians Ha.b so that we can take the limit ipuisc — )■ 0, /2 — >■ oo with I2ipuise = 
7r/2. After integration of the Heisenberg equations of motion during tpuise, it is found 
that the fields arc transformed by the 7r/2 pulse as follows: 

V;a(0+) = -^[^.(0-)-V'6(0-)] (9) 

V>b(o+) = -i=[^,(o-) + V'b(o-)] (10) 

We are interested in the squeezing and quantum correlations that develop during the 
non-equilibrium dynamics following the pulse for i > 0. 



2.2 Our Approach 

The problem of the scaling of the squeezing for TV — >■ oo in the multimodc case 
implies the solution of the non-equilibrium quantum dynamics for a large number 
of atoms and a large number of modes. We cannot solve exactly this problem even 
numerically. However, what can be solved exactly on a computer is the "classical 
field equivalent" of our problem. We then adopt the strategy summarized in Table [TJ 
We use the classical field model to (i) perform numerical experiments and (ii) test 
a perturbative solution that we can generalize to the quantum case. The quantum 
perturbative solution is then used to get quantitative predictions on the real physical 
system. 



^ The exact relation between the bare coupling constant and the effective coupling constant 
is 1/5 = 1/30 + /pBz where FBZ=h^/Z, 
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Table 1. We cannot solve exactly the problem in the quantum case but we find an analytical 
perturbative solution. We check our perturbative approach in the classical case where we 
can solve the model exactly (numerically). 



Quantum field 
model 


solution 
available ? 


Classical field 
model 


solution 
available ? 




no 


H[Mr),Mr)] 


yes 

t 

yes 


Perturbative 


yes 


Perturbative 



2.3 Classical field model 

The classical field model |21I22| is obtained by replacing the quantum fields with 
classical fields in the Hamiltonian 

"01: V'a, V'l, V'h ^p*,1pa,^t''^b- (H) 

In the equations of motion the commutators are then replaced by Poisson brackets. 
The classical field model is useful when the interesting physics is given by low-energy 
highly populated modes |23I24I25] . For our classical field simulations, we assume that 
this is the case in the equilibrium state before the pulse, with all the particles in state 
a. The initial field ipa'^ then randomly samples the thermal equilibrium classical field 
distribution for the canonical ensemble at temperature T 

Pel = |e-'^^ (12) 

where H is the classical Hamiltonian, which is a discrete version of the Gross- 
Pitaevskii energy functional: 

^ = E ^^^^^ + ^dvY,rair)ra{r)Mr)Mr) (13) 

k r 

For the initially empty state &, inspired by the Wigner quasi-probability distribution 
of the quantum density operator and the truncated Wigner approach j26l27l28l29] 
we represent the vacuum by a classical field "0^"' having in each mode independent 
Gaussian complex fiuctuations of zero mean and variance 1/2: More precisely, we set 

^k*"* = X + iY where the independ 
Gaussian probability distribution 



^k*"* = X + iY where the independent real random variables X and Y have the same 



At t = the fields are mixed by the pulse according to ©-(ITU]). At later times, the 
fields, il)a and '06 evolve independently according to the discrete non-linear Schrodingcr 
equation {v = a, b) 



(15) 



^ In this classical limit, the substitution go — >■ g is required, since the difference between 
l/go and 1/g is due to quantum fiuctuations. 



6 



Will be inserted by the editor 



where the discrete Laplacian A has the plane waves exp(ik • r) on the lattice as eigen- 
vectors of eigenvalues —k^. The lattice model automatically provides a momentum 
cut-off to the classical field model, corresponding to the boundaries of the first Bril- 
louin zone FBZ=:[— tt//, 7r//[3. The various observables have a more or less pronounced 
dependence on the cut-off. Here, guided by our analytical results (see section S]) we 
choose the cut-off such that, in the thermodynamic limit, the non-condensed density 
for an ideal gas in continuous space in the Bose condensed regime (zero chemical 
potential) is exactly reproduced by the classical field model: 

f d^k 1 f _f^kBT 

U (2^)3 ePE'^ - 1 ipBz (2^)3 Ek ^ 

where Ek = fi^k^ /2m is the kinetic energy, the mode occupation numbers are given 
by the Bose formula for the quantum case and by the equipartition formula for the 
classical case. As revealed by the change of integration variable K = Adsk, the con- 
dition ([TS]) is an equation for l/XdB with AdB = [27rft^/(mfcBr)]^/^ the thermal de 
Broglie wavelength. The integrals on both sides of (flB)) can be calculated analytically, 
see e.g. [30] for the integral in the right-hand side, and the usual factor C(3/2) appears 
in the left-hand side (where C is the Riemann Zeta function) . The condition then 
gives E™'^^ ~ 2.695kBT with E™'^^ is the maximal kinetic energy on the grid, here 



3 Numerical Experiments 
3.1 Dimensional analysis 

As specified in subsection l2.11 when the lattice model approaches the continuous space 
physics for our observable (the spin squeezing), the physical parameters of the model 
are the atom mass m, the effective coupling constant g characterizing low energy 
binary interactions between the atoms, the temperature T, the total atom number N 
and system volume V: 

-, g, kBT, N, V (17) 
m 

The spin squeezing parameter optimized over time, ■Cmim ^ dimensionless quantity. 
It is therefore a function of the independent dimensionless combinations that we can 
form from the ensemble (fT7|l : 

^l-^^fiN,^,^) (18) 

Here \/ pa? is the "small parameter" such that ^ pa^ ^ 1 characterizes the weakly 
interacting limit, and pg is the mean field chemical potential of the gas (at T = 0). 
The same dimensional analysis and the same general form of ^^jj^ hold for the classical 
field model. 



* In the classical field simulations, to ensure maximal ergodicity, we did not use a cubic 
lattice. We used the following aspect ratios for the quantization box, : Ly : = y/2 : 
{1 + V^)/2 : ^/3, and we discretized the positions with the same even number of points Wmax 
along each direction. This corresponds in condition (|16|l to a slightly asymmetric FBZ. 




10"* 10^ lo'' 10^ 

N 



Fig. 1. (Color online). Minimal squeezing parameter (that is, minimized over time) 
for four different temperatures and increasing system sizes in the thermodynamic limit 
TV — >■ oo, p,g,T = constant. Squares: classical field simulations result. ksT / pg = 
1.13(a), 0.78(6), 0.50(c), 0.28(d). For all the points i/p^ = 1.32 x 10"^. The horizontal 
dashed lines are analytical results in the thermodynamic and weakly interacting limit H21|) 
that are the classical field equivalent of (|62|) . The grid sizes (number of points per direction) 
are n^ax = 12, 16, 20, 32, 36 for (a), n,^ax = 10, 12, 16, 32, 40 for (b), n^ax = 12, 16, 24, 36, 40 
for (c) and n^ax = 6, 8, 12, 16, 24, 32 for (d). 



3.2 Existence of a thermodynamic limit for 

We have performed classical field simulations, increasing the system size in the ther- 
modynamic limit 

^ oo ; V ^ QO] p, g, T ^ constant , (19) 

In Fig|T]we show the result for four different temperatures. The squeezing parameter, 
minimized over time, converges to a finite value. According to the general form (jl8p , 
^min then depends on ksT/pg and yfpa^. The first parameter ksT/pg is varied in 
FiglU whereas the parameter y/ pa^ defining the weakly interacting regime is main- 
tained constant in that figure. Note that for curves (a)-(c) the limit is already almost 
reached for = 3 x 10^ while a larger system is needed for the lowest temperature 
curve (d). 



3.3 Weakly interacting limit of 

Starting from a point that is already at the thermodynamic limit for each temperature, 
we have performed other simulations, going more deeply in the weakly interacting 
limit p oo, g —i' 0, pg,T = constant [JT] 0. In this limit the small parameter \fp^ 
tends to zero. In Figl^l for a fixed value of ksT / pg, we show the squeezing parameter 
divided by \fp(^ as a function of a rescaled time, for various values of \/ po? . It is 



^ In our simulations (with finite size systems) we increase A*' and decrease g while V , T 
and Ng are fixed. 
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Fig. 2. (Color online). Squeezing parameter divided by \ 
time, in the weakly interacting limit: p — s> oo, g — > with pg,T — constant. Solid lines: 
classical field simulations. Parameters: ksT / pg = 0.5 and Jimax = 32 for all the curves, (a) 



iV = 3 X 10* 



= 1.32 X 10"^; (b) N = 10" 



= 3.96 X 10"''; (c) iV = 3 X lO" 



■ = 1.32 X 10"^ (d) iV = 6 X 10^ 



— 6.59 X 10 . Horizontal line: analytical result 



apparent that the rescaled minimal squeezing is nearly constant in the figure. For 
weak interactions, ^'^^^/ \/ po? is thus a function of UbT / pg only. 



3.4 Spin squeezing as a function of the temperature 

From the numerical experiments described in 13.21 and 13.31 we conclude that, in the 
double thermodynamic plus weakly interacting limit, the squeezing parameter mini- 
mized over time is equal to \l~pa^ times a universal function of ksT / pg: 

^^:„-V^^(^) (20) 

We show the universal scaling of the minimal squeezing parameter in the thermody- 
namic limit in FiglSl where we collect the results of several simulations for different 
temperatures and interaction strengths. In the following section we shall develop the 
analytical theory that gives the explicit expression of the function F appearing in 
((20| . The analytical result for the classical field theory is represented as a full line in 
FigH It reads: 

Here p = N/V is the total density, Sk = Uk + Vk, s[°^ = t/f.^' + V^'"^ are Bogoliubov 
functions defined in ((30)) and ([28| . n]^^ = ksT/e]^^ are the equilibrium occupation 
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kgT/pg 

Fig. 3. Minimal squeezing parameter ^^j^ divided by yfpcfi as a function of kaT / pg. Sym- 
bols: classical field simulations with \J pa? = 1.32 x 10~^ (squares), 1.94 x 10~^ (disks), 
and y/'pefi = 4.17 x 10~* (triangles). Solid line: Analytical classical field result (|2ip . Dashed 
line: Quantum result ((62]). The grid sizes, for increasing ksT/pg are rimax = 32, 36, 40, 36 for 
squares, rimax = 24, 32, 32 for disks and nmax = 32 for triangles. The initial non condensed 
fractions in component a before the pulse are {Nnc)/N — 0.01,0.02,0.04,0.07 for squares, 
{Nnc) /N = 0.01, 0.02, 0.04 for disks and (iVnc) /N = 0.01, 0.02, 0.05 for triangles. 



numbers of Bogoliubov modes before the pulse with e^°^ = [Ek{Ek + 2pg)Y/'^ and the 
integral is restricted to the first Brillouin zone FBZ with G [—k™^^,k^^'^^[ for the 
three directions v = x,y,z. We note a very good agreement between the analytics 
and the simulations for all the points. For comparison, we also show the analytical 
result for the quantum field (j62p . in the zero-lattice spacing limit, as a dashed line. 
Note how well the classical results with the cut-off prescription (|T6)l reproduce the 
quantum analytical results for ksT > pg. 



4 Analytical results 

We want here to address the fully quantum case and calculate analytically the func- 
tion F{kBT / pg) = \/ pd^ of (^0]) using Bogoliubov theory. In 14.11 we present 
a modulus-phase reformulation of the Bogoliubov theory generalized for a bimodal 
condensate; in 14.21 and 14.31 we derive the expansion of the spin squeezing parameter 
in the thermodynamic limit, and we discuss the final analytical results in 14.41 

4.1 Modulus-phase Bogoliubov formalism for bimodal condensates 

We generalize here to two-components the J7(l)-symmetry preserving Bogoliubov 
theory of [21], see also [J. As spin squeezing in bimodal condensates is due to phase 
dynamics, we rephrase the theory in terms of the relative phase operator of the a and 
h condensates. This is a crucial step that allows us to obtain results by a pcrturbative 
expansion in powers of that relative phase. The relative phase is simply the difference 
of the condensate phases 0a, h introduced as hermitian operators conjugate to the 
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condensate atom numbers |32l33j . This is a valid description in the subspace excluding 
the vacuum state (no particles) for the condensate modes, which is sufficient here: as 
we suppose initially T ^ Tc, both condensates are highly populated after the mixing 
pulse with (No) = (Ni,) = N/2. Some useful expressions involving the Bogoliubov 
amplitudes and some commutation relations are given in Appendix 1X1 

Within each atomic internal state a or b, we split the bosonic field operator into 
the condensate and non condensed modes contributions: 

ao ' bo - 

'<Pa = + '0a_L , tPb = + ^6_L (22) 



For the annihilation operators in the condensate modes we introduce the modulus- 
phase representation 

ao = e'^"\/^ ' \^a^.k\=i (23) 
foo^e^'^V^ ' [^bo,^b]-i, (24) 
while for the non condensed modes we introduce the number conserving fields [31134] 

= e-'^"7^,i ; i6=e-'^'^b± (25) 

that we expand over Bogoliubov modes with amplitudes Cak and Cf,k- To perform this 
expansion, one has to distinguish the time before the pulse and the time after the 
pulse. 

Before the pulse, at time t = 0~ , all the N atoms are in the internal state a, with 
a zero-temperature mean field chemical potential ^ = pg, and we expand the number 
conserving field Aa{t = 0^) over the Bogoliubov modes: 



^i°^-EN°'4l^ + ^i°^(^iV]^ (26) 



Here the exponent over the operators indicates that the operators are considered 
at time t = 0~ . The bosonic operator cf^ annihilates a Bogoliubov quasi-particle of 
wave vector k, and eigenenergy 

ei"^ ^[E,{E, + 2pg)]^/' (27) 

where we recall that Ek ~ f?l? j (2m) is the kinetic energy contribution. The cor- 
responding amplitudes of the Bogoliubov mode, correctly normalized as — 
= 1, are given by 

1/4 



s 



f .t/r + ^r-rrTHT^-flT^) (28) 



Before the pulse, the field in the b state is in the vacuum state, so the modal expansion 
is performed over the usual single particle plane wave eigenbasis. 

After the pulse, at < > 0+ , the particles are on average equally distributed among 
the two internal states a and b, (Na) = {Nb) = N/2. There are now two condensates, 
with interaction constants gaa ~ gtb ~ 9 among atoms of same internal state, and no 
interaction [gab = 0) among atoms of different internal states. Similarly to p6p. we 
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perform after the pulse the modal decomposition of the number conserving fields in 
each internal state a ~ a,b: 

Aa = J2 + Vkcl^u] ^ (29) 



where c^k annihilates a Bogoliubov quasi-particle of wave vector k in internal state 
a. The Bogoliubov mode amplitudes Uk, Vk and the eigenenergies do not depend 
on the internal state and are deduced from the ones at t = 0~ by replacing the mean 
field term pg by pg/2: 

1 / E, 



Sk^U, + Vk = = — ^ (30) 

Uk-Vk \Ek+pgJ 

ek^[Ek{Ek+pg)Y'\ (31) 

Note that this involves an approximation: In principle, the Bogoliubov modes in in- 
ternal state a = a or 6 depend on the actual number of particles N„ in that state, 
which has small ~ 1/VN relative fluctuations since, after the pulse, N^j has a bi- 
nomial distribution peaked around 7V/2. Taking into account this effect changes the 
mathematical structure of the theory, since the Bogoliubov amplitudes Uk and Vk, 
and thus the quasi-particle annihilation operators Co-k, would then depend on the to- 
tal number operator Ncy in state a, which is beyond the scope of the present work. We 
nevertheless verified numerically on the complete Bogoliubov theory (in the classical 
field model) that this fixed-Bogoliubov-mode approximation is extremely accurate 
both at short and long times, introducing (for the typical parameters considered in 
our figures) a relative error on lower than 10^^ comparable to our statistical error 
bars with 10^ realizations. Another important point is that the Bogoliubov quasi- 
particles are not at thermal equilibrium after the pulse, so that, for example, their 
mean occupation numbers are not given by the Bose formula. At the level of the Bo- 
goliubov approximation, the quasi-particles do not interact and cannot thcrmalize, 
the corresponding quasi-particle creation operators evolve in Heisenberg picture with 
simple phase factors: 

c.k(t) = e— ''*/''c,k(0+) (32) 

for a = a,h. The validity of this no-thermalization approximation is discussed in 
subsection 14.41 

To express the evolution of the particle annihilation operators ao and feo in the 
condensate modes, we use the modulus-phase representation (|23l [^^ . For the mod- 
ulus, one simply uses the conservation of the total atom number in each internal 
state, 

Nao ^N,-N,^^N,-Y, dVAlA, (33) 

r 

SO that iVcro can be expressed in terms of the quasi-particle operators Cak- For the 
phase, we use within each internal state a ~ a,b the equation of motion derived for a 
single component in }16) and truncated at the level of the Bogoliubov approximation: 

j^0„ --x(n.-1)-^Y1 dV {Al + 2it + it2) (34) 



where we have introduced x = 5/(fi-^)- Replacing the operators A„ by their modal 
expansion (j29p . one gets contributions that do not oscillate in time, and contributions 
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such as C(jkCcr-k that oscillate in time (at the frequency 2ek/h in the example). As 
we shall need the value of the phase operators at long times t (typically efci/fi- ^ 1) 
rather than the value of its derivative, we argue as in jl6| that the oscillating terms 
in p4p . after temporal integration, give a negligible contribution to the squeezing 
parameter. We have checked this approximation analytically: The expression for 
fully including the oscillating terms in the phase difference operator is given in the 
Appendix |E1 it corrects the approximate expression (|60p of by typically a sub- 
percent effect at intermediate times, and by a vanishing amount at large times. These 
oscillating terms in 9a — are thus of little physical relevance for the spin squeezing. 
Keeping only the non-oscillating terms gives for the relative phase operator of the 
two condensates at time t: 

{Oa - k){t) = [Oa - Ob){Q^) - Xt [Na - iVf, + (35) 

'D^Y.^Uk + Vuf{nai.~n^].){Q+) (36) 

k^^O 

where we have introduced the quasi-particlc number operators no-k = c^^j^Co-k, which 
are constants of motion in the Bogoliubov approximation, see ([32]). The multimode 
contribution T) to the relative phase ([55]) will play a central role in what follows. 
It is indeed because of this term (neglected in the usual two-mode models) that the 
squeezing parameter is bounded from below by a non-zero value in the thermodynamic 
limit. 

The last step is to relate the various operators at time t = 0+ to their values just 
before the pulse. Since the state of the system is known at i = 0~, this fully specifies 
the "initial" conditions for the time evolution of the operators after the pulse. The 
derivation and the more precise results are given in the Appendix |B] Here we give 
the main conclusions. The initial value of the condensate phase difference is, in the 
large N limit: 



>)r) = ^{(»sr"'('i°'-'-)}+«(i^'^ 



(37) 



where b^Q^ = e~^^"'^bQ^ and { , } stands for the anticommutator. This results from 
the coherent mixing of the initial condensate amplitude with the vacuum noise fluc- 
tuations in the initially empty internal state 6, in the same spatial mode as the initial 
condensate in state a, that is in the plane wave with zero wave vector. Similarly, the 
quasi-particle annihilation operators just after the pulse are coherent superpositions 
of the initial field fiuctuations, mainly thermal fiuctuations A in a and only vacuum 
fiuctuations B in b. In the large N limit, 

c,.(0*) = ^ + o(-^) (38) 

where the — sign is for cr = a and the -I- sign is for a — b. The expression of in 
terms of t = 0~ operators of the a internal state, and the expression of J5k in terms 
of i = 0~ operators of the b internal state, naturally appear in the calculations of 
Appendix |B] 

ik ^ (f/.C/f - V.V!^') cit^ + (C/.Ff ) - y.C/r ) ci°il (39) 

i3kE.C/.e-«'"'6f -^.^^'""^LT (40) 

In Appendix[Cl it is pointed out that these number conserving operators obey bosonic 
commutation relations and all their second moments are explicitly evaluated. 
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4.2 Double expansion method in the thermodynamic and weakly interacting limit 



We shall now apply the Bogoliubov theory developed in section I4?l1 to the calculation 
of the time-dependent squeezing parameter defined in ([3]). For the configuration 
that we consider, symmetric under the exchange of a and 6, the mean spin is always 
aligned along x. The minimum transverse spin variance is then 



{sD + {sD ^{{si)-{sw + {{s..SyW 



(41) 



where { , } is the anticommutator. From the definition ^ it appears that 5*2 is a 
constant of motion, its variance can thus be evaluated just after the pulse, at i = 0+. 
According to (|9ll0p . the pulse applies to the collective spin a rotation of angle 7r/2 
around y axis, so that 

5.(t>0) = -^(°) and {SD^^ (42) 

To obtain the spin variance {Sy) and the spin correlation {{Sy, Sz}), the challenge is 
to determine, as a function of time, the operator Sy, or equivalently the antihermitian 
part of the operator S+ introduced in ([T|) (in its discrete version for the lattice model), 
since Sy = {S+ — h.c.)/(2i). In the expression of 5+, one applies the splitting (|22l) of 
the bosonic fields in the condensate and the non-condensed contributions, one uses 
the modulus-phase representation (j23l for the condensate part and one introduces 
the number-conserving fields for the non-condensed part, to obtain: 

S+ = e-'(^»-^') + p] with (43) 



F = ^/i^^^^+l)^o - Y + E dVAlAh (44) 

r 

Guided by the numerical experiments in section |21 we have developed a systematic 
double expansion technique to determine (Sy) and {{Sy,Sz}). The two small param- 
eters controlling the large system size limit [i.e. the thermodynamic limit (jl9p ] and 
the Bogoliubov limit are 

Csize = and CBog = — ^ (45) 

where (^„c) = {N^f/) is the mean number of non-condensed particles in the initial 
state, which is indeed much smaller than for a weakly interacting gas at T ^ Tc. 
For the Bogoliubov expansion, we will keep terms up to order one included in eBog] 
keeping higher order terms would not be consistent with the use of the quadratic 
Bogoliubov Hamiltonian. To determine the required order of the large system size 
expansion, we note that (Sx) in the denominator of ^ remains close to its f = 0^" 
value N/2 over the relevant time scales (that are finite in the thermodynamic limit 
with pgt/h ^ N^^^), so that « 4:ASj_ ^^^/N. To have a vanishingly small error on 

in the thermodynamic limit, we will keep in {Sy)/N and {{Sy, Sz}) /N terms up 
to order zero included in esizc, that is we can neglect the contributions that tend to 
zero when egize — >■ 0. 

The systematic technique to determine the order of an operator is to estimate its 
mean value and its standard deviation in the quantum state of the system. This is 



14 



Will be inserted by the editor 



safer than a simple guess, in particular when the operator has a vanishing expeetation 
value. A relevant example is the operator T) defined in p6p . that will play a crucial 
role in the best achievable squeezing. After a superficial look at p6p . one may believe 
that T) scales as N in the thermodynamic limit, since it involves a sum over all modes, 
and that it scales as eBog in the Bogoliubov limit since it involves the quasi-particle 
number operators fi„\^. However, it is actually the differences n^k — "-bk that matter, 
and for the particular state resulting from a 7r/2 pulse applied on a gas initially in 
the a internal state. As a consequence, the expectation value of T) is zero, it is the 
variance of V which scales as iV, so V scales as N'^^'^ in the thermodynamic limit. 
To determine its scaling with CBog, we keep the leading term (|38D in the value of the 
quasi-particle annihilation operator just after the pulse, to obtain 

P~ -^(C/fc + H)2(4Bk + S^ik) (46) 
This scales with N as N^/"^ since each term has a zero mean. Intuitively, this scales 

1/2 " . 

with eBog as CBog- is of order unity since it corresponds to vacuum field fluctuations 

" . 1/2 

in the initial empty state 5, and is of order e^^^ since it corresponds to the initial 
non-condensed field fluctuations in state a. We thus conclude that 

^^(iVeBog)^/' (47) 

This is confirmed by the correlation functions of A and B given in the Appendix [Cl 
that allow an explicit calculation of {'D'^)/N, which is indeed « CBog, see Eas. (j61l62p . 

The same analysis is applied to the various operators appearing in the antihermi- 
tian part of (|43p . writing for simplicity F = Fn + iFj, where Fji and Fj are hermitian. 
It is found in the Appendix [D] that 

^^-^^^J^ (48) 
Fi « (iVeBog)'/' (49) 
Fb. ~ iVcBog ± (A^EBog)'/' (50) 

Contrarily to the first two operators, Ffi has a non-zero expectation value, and the 
writing of its estimate in ([50)) corresponds to its mean value plus or minus its standard 
deviation. It turns out that the fluctuations of Fr are negligible so that the operator 
can be replaced by its mean value. The weak value of the phase difference operator 
shows that its exponential in can be expanded to first order, which substantially 
simplifies the calculations. The final result, up to zeroth order included in egizc and 
up to first order included in eBog, is 

h2%^) (51) 

{{Oa - e^f) ( J + {F„)^ \{{FiA - k}) (52) 
- + {Fr) \ {{Na - m, Ba - Bb\) (53) 





4(. 






N 


N 




1 


N 


^ ^2N 



Note that these expressions hold independently of the approximation performed in 
the phase difference operator (that neglects the oscillating terms). 
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4.3 Results of the expansion method for ^^(i) 

To obtain an expression for in the double thermodynamic and Bogoliubov limit, 
it remains to explicitly evaluate (j51l52l53|) and to insert the result in pi4ip . The 
required operators have their expression given in useful form in the Appendix |D1 and 
by P5I37I46P for the phase difference operator. 

We have first evaluated (|51l52l53p at time t = 0+, and we have checked that one 
recovers the exact relations (S'^)(0+) = N/2, (S'2)(0+) = 7V/4 and {{Sy, S^}){0+) = 0. 
At finite time, squaring ([35]) . one realizes that the crossed term, linear in time, has a 
zero expectation value since 

{{0a - Ob){Q+){Na - Nb)) = {{Oa - Ob){0+)'D) = (54) 

This is due to the fact that the phase difference operator at t = 0+ is proportional to 
the antihermitian part iy of e~*^" '^0*^ whereas Na — Nt only involves the hermitian 
part X and T> does not depend on that operator. As a consequence, the spreading 
of the relative phase is purely ballistic within our Bogoliubov approximation, that it 
ignores the interactions among the quasi-particles 

{{Oa - e,f) = {{0a - ^fc)'>(0+) + {xtf{{Na - N, + Vf) (55) 

(the inclusion of these interactions within a quantum kinetic framework introduces a 
diffusive component in the phase spreading jl8j). Another simplification takes place, 
for similar reasons, 

{Fi{Na- Nb+T))) (^i^ (56) 

so the last term of ([5^ reduces to its < = 0+ value. Using l^^ . and the fact that 
(Fr), (I?^) and {{Sz.V}) are all « iVesog, and neglecting contributions in e|og, we 
finally obtain 

{^1) - {SD _ ^2 ( 1 , \^ I \ , ^g^^ 

(58) 




N 

where we have introduced a dimensionless "time" 

-w-f ('+ '™"';"'-'^" ) 

that is slightly renormalized by a time dependent contribution of order esog- The 
expectation values {Fji{t)) and {{Sz, 2?}) appearing in (j59p are given in the Appendix 
IdI they are respectively uniformly bounded in time and time independent. Note that 
the long-time behaviors of (j57l58p were expected: As the squeezing dynamic occurs, 
Sy indeed grows quadratically in time while Sz stays constant. 

Finally, expanding (|4ip up to order one included in esog at any fixed r gives the 
squeezing parameter as a function of time in the thermodynamic limit (in particular 
for T < iVi/2): 

^2 ^ ^ N ^ ^ ^ 1_L^ (60) 
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Fig. 4. (Color online) Black solid line: Squeezing parameter in the thermodynamic and 
weakly interacting limit (|60|) as a function of the renormalized time r (|59|l . Parameters: 
ksT/pg — 1 and \fp(^ ~ 10^'^. Black dashed line: asymptotic value (|6ip . Red solid line: 
two-mode model prediction. 



As we show in FigUl the squeezing parameter ^^(i) decreases in time essentially as 
in the two-mode model [this is the first term in the right hand side of (pO)) without 
the {Fii)/N term, see e.g. equation (52) in [TT])] until a renormalized time r w 
l/^EBog ^ 1 is reached, when the multimode effects [the second term] start limiting 

the squeezing. At such times, the contribution of {T)'^)/N dominates over the one of 
{Fii)/N by a factor « l/eg^g, which shows that {'D'^)/N constitutes the real actor in 
the process of squeezing limitation. 



4.4 Minimal squeezing and best squeezing time 

From the central result (j60p of the previous subsection, it appears that the minimal 
squeezing ^^^^j^^ is reached in the thermodynamic and Bogoliubov limits at "infinite" 
time and it is given by 



An explicit calculation gives: 



Sn 



N 



(61) 




(62) 



where we have introduced the mean occupation numbers nf^ = l/[exp(/3e)."'') — 1] 
of the Bogoliubov quasi-particles in the initial (t = 0~) thermal equilibrium gas in 
internal state a. The prediction of (p^ is shown as a dashed line in FigE] and as a 
full line in FiglS] In Figl5]wc show that the minimal squeezing parameter given by 
([5^ is always lower than the non condensed fraction {N-nc)/N where (iVnc) = (-^ix) 
is the mean number of non condensed atoms in component a before the pulse: 



(0) 



(0) 



,(0) 



(0) 



(63) 



Asymptotically, for ksT ^ pg (but always T <^ Tc) the minimal squeezing parameter 
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Fig. 5. (Color online). Minimal squeezing (black solid line) ^^i^ and non condensed fraction 
(red dashed line) {Nnc)/N, both divided by \/pa^, as a function of kBT/pg. 



^^ijj reaches the non condensed fraction. In the opposite limit, for ksT/pg — > 0, the 
squeezing tends to a constant value. 



f 2 (T=0) 



fV2-^ln(V2 + l) 



0.02344 (64) 



Although non zero, is very small for practical purposes. Indeed pa^ < 10 ^ 

in present squeezing experiments so that predicts Cm-n""'' ^ 2 x 10~^. 

The fact that the minimal squeezing is obtained at infinite time is a limitation of 
our Bogoliubov approach, that neglects the interactions between the quasi-particles 
and effectively assumes that the corresponding collision time is diverging, see discus- 
sion below. However, since the numerical squeezing curve C^(i) is quite flat around 
its minimum (see Fig [5]), it suffices in practice to determine a "close to best" squeez- 
ing time tjj defined as ^^(t,,) = (1 -I- ?y)Cmini where 77 > 0. Then, according to 
expanded for large r up to order included, t.,j is finite and given for <C 1 by 

The "close to best" squeezing time for 77 = 0.2 is shown in FiglHlas a full line 

and compared to simulations (filled symbols). 

An important issue is that of thermalization that brings the system back to equi- 
librium after the pulse. Thermalization is neglected in the Bogoliubov theory and in 
our analytics but it is included in the classical field simulations. It is thus possible 
to reach ^-^ (1 + ?/)^^jjjj only if given by ((65|) is shorter than the thermalization 
time: 

trj < ithorm (66) 

We show the thermalization times, that we extract from the classical simulations as 
explained in [S], as empty symbols in FiglHl For the points we considered they are 
indeed longer than the close to best squeezing times and the condition (|66)) is satisfied. 
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Fig. 6. Universal scaling of the "close to best" squeezing time trj with -q = 0.2. Filled 
symbols: classical field simulations with a/ pa^ = 1.32 x 10~^ (squares), 1.94 x 10^'^ (disks) 
and ^/ pa^ = 4.317 x 10~* (triangles). Solid line: analytical result (|65p using (|2ip . The empty 
symbols (squares, circles and triangles) show the thermalization times ftherm extracted from 
the simulations as explained in [S]. Parameters are as in Fig[31 



We can also estimate ttherm from Landau-Beliaev damping rates of Bogoliubov modes 
|17l35j . The damping rate of the mode q is 

"Shoal 

where the healing length ^heai is defined by pg = h-^ / (2m^[^^^j) and the rescaled Landau 
and Beliaev damping rates and are dimensionless functions of ksT / pg only, 
given e.g. in equations (A7) and (A13) of [T7]. Concentrating on the pre- factor in 
(|57)l . for fixed ksT/pg, this gives 



pgtthcim ^ pg 
h hFq 




(68) 



The scaling with (pa^) of the thermalization time is shown in FigEl On the other 
hand, the close to best squeezing time i,, scales as 



ft (pa3)^/^ tth 
and (j66p is satisfied in the weakly interacting limit 



so that ^^(x(pa3)i/4 (69) 



crm 



5 Physical interpretation 
5.1 Limit to the squeezing 

In the spin squeezing scheme with condensates that we consider, the useful quantum 
correlations are built through mean field interactions that introduces a phase shift 
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Fig. 7. Scaling with (pa'*) ^^'^ of the thermahzation times extracted from the classical field 



simulations, 
(triangles). 



= 1.32 X 10"^ (squares), 1.94 x 10"^ (circles) and ^Jpcfi = 4.317 x 10" 



for each atom that depends on the collective variable Na — N),. In the collective spin 
picture, we can say that in a given realization of the experiment, the component Sy 
becomes an enlarged copy of Sz so that correlations build up in the Sy-Sz plane, 
orthogonally to mean spin. To explain this fact in a simple reasoning, we can identify 
Sy with the condensate relative phase: Sy ~ ^(^a — Ob) and look at equation (|35|) 
that we rewrite here replacing x by its expression x = 9/{KV): 



Na-Nb+V 



(70) 



Initially, at t = 0, the phase difference (6*0 — 6'f,)(0+) is of order l/vN and Na — Nb is of 
order As soon as pgt/h ^ 1, the time dependent term in ([701 dominates over the 
initial condition. In the absence of the multimodc contribution to the phase difference 
(i.e. for V = 0), 9a — 0b and thus Sy become an enlarged copy of Sz = {Na — Nb)/2. 
This is the scenario in the two-mode theory. Correlations between Sy and Sz becomes 
perfect in the long time limit. In this case there is no limit to the squeezing and 
Cmin ^ when N ^ 00. On the other hand, in the presence of T> this is not possible. 
Looking at squeezing in the long time limit where \Sy \ 3> \Sz\ and keeping only the 
leading (O"^) order in 1/Mn equation (pO)) . we can write 



{s'y){s'z)-{{Sy,Sz}/2y 



{S'){S^z 



N 



(71) 



The only contribution left in ^^j^^ is the variance of V that is the part of 6a — Ot that 
is not proportional to Na — Nb- But what is the physical origin of T), given by ([55]) ? 
It comes from the fact that the mean field interaction for a condensed atom with and 
another condensed atom or with an atom in an excited mode is not the same. This 
is particularly clear in the Hartree-Fock limit where 14 — >■ and Uk — >■ 1- In this case 
2? reduces to Nai_ — Nb± that is the non-condensed atom number difference. In this 
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limit we have 



HF 



hN 



(72) 



the factor 2 is the Hartree-Fock factor that doubles the effective strength of the 
condcnsate-non condensate interaction with respect to the condensate-condensate 
interaction. 



5.2 Squeezing of the condensate mode 

In FiglSl for two temperatures: fcsT ^ pg and ksT < pg, we compare the squeezing 
of the total field that we have been considering so far, with the squeezing of the 
condensate mode constructed with a spin operator involving the condensate mode 
only: Sqx + iSoy = aJ,6o and 5*02 = Nao - Nw- 

This is the squeezing that would be obtained by "selecting" only the condensed par- 
ticles for the squeezing measurement. Besides the classical field simulations, in Figl8^ 
and FiglSb we also show as dashed curves results obtained in the Bogoliubov approx- 
imation Q Clearly <C in both graphs. Particularly striking is the case in Fig|S)D 
where Comin/^min — ^0 whilc the non condensed fraction is only {Nnc)/N = 0.02. We 
explain here why this is the case. In order to have condensate squeezing we need to 
build up correlations between Soy, that is still proportional to 9a — Ot, and Sqz- Ac- 
cording to ([7D| , at long times 9a — 9b differs from Sqz by the quantity Na± — Nb± + V 
that prevents the correlations to become perfect at long times. Indeed, we find at long 
times and for a large system that 



Var 



(Nal^ - NbA) + 



V 



m) Tr (74) 



N 

The evaluation of the minimal achievable values of from (j74p is more involved than 
for because the quantity Nai. — Nb±, contrarily to V, is not a constant of motion. A 
detailed discussion is beyond the scope of this paper, but one can give simple reasons 
explaing why the minimal value of is numerically found to be much larger than the 
one of Let us first forget about the time dependence of Na± — Nb± and evaluate 
it at time t = 0+. It is found that the variance of Na± — Nb± at t = 0+ is simply 
{nI°1), that is the mean number (A^nc) of non-condensed particles before the pulse. 
In the Hartrcc-Fock limit ksT ^ pg, V reduces to Na± — Nb±. One then expects 
that the minimal is four times the non-condensed fraction, that is four times larger 
than ^^i,!- In the low temperature regime fc^T <C pg, (V^) <C {Nnc), so the ratio of 
condensate to total field squeezing is expected to be even larger. 

Let us now take into account the time dependence of Na± — Nb±. This makes the 
situation even worse for the condensate mode squeezing: Whereas the operator V has 



^ Before the pulse, we start with a thermalized field sampling (|12|l . After the pulse, 
we evolve the condensate phase with the classical equivalent of H34I) . also performing in 
that equation the approximation of neglecting the oscillating terms, and the Bogoliubov 
amplitudes with (|32p . The condensate atom numbers are obtained by the classical equivalent 
of (O. 
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normal fluctuations with a variance scaling as the volume V of the system, it is found 
that Nai_ — Nh± has anomalous fluctuations at long times. If one replaces discrete 
sums over k by integrals in the expression (|117p of the variance of Na± — Nii±, as 
usual in the thermodynamic limit, one finds converging integrals but the variance 
diverges linearly in the long time limit: 



Var(jV„^- jVb^) 
TV 



therm. lim. 



This result obtained within the Bogoliubov theory of course fails for times larger 
than the thermalisation time. In practice, it is more physical to consider a finite 
size system. One then finds that, in the long time limit, the variance is dominated 
by the contribution of the low-fc terms in the sum over k in (|117p : The variance 

of Na± — Nii± grows from its extensive t = value {N^^^) to a super-extensive 

value ~ 10{kBT/pg){mcL/2Trh)'^U 

in a time scaling as i/c where L = V^^^ is the 
system size and c = ^Jpglm is the initial sound velocity. After this time, the variance 
oscillates around the super-extensive value with a period again scaling as L/c, given 
by the fundamental Bogoliubov modes in the box. This explains the oscillations of 
and its temporal mean value in FiglB}). In FiglB^ the oscillations are less visible as 

the anomalous contribution to the variance is smaller than its i = 0"*" value (iV^*^^). 



6 Conclusion 

We have shown that, in a multimode theory, the spin squeezing that can be obtained 
dynamically using interactions in condensates is finite in the thermodynamic limit. 
This is contrary to the results of the currently used two-mode theory that predicts 
an infinite metrology gain for iV — >■ oo. Using a convenient reformulation of the 
Bogoliubov theory, we could calculate the temperature and interactions dependent 
limit of the spin squeezing parameter analytically for a spatially homogeneous system. 
We performed non perturbative classical field simulations to test our analytical results 
including interactions among Bogoliubov modes and thermalization that are neglected 
in the perturbative treatment. 

At temperatures /cbT <C pg the limit that we find for the squeezing parameter 
optimized over time, C^jij,, is very small and in particular much smaller than what is 
currently measured in present experiments. Nevertheless it represent the fundamental 
limit of this squeezing scheme and we hope that the temperature dependent limita- 
tion to spin squeezing will be soon within reach of experiments. We explained that 
the physical origin to this limit of the squeezing lies in the difference of mean field 
interactions between condensate-condensate and condensate-non condensed atoms. 

A.S. acknowledges useful discussions with J. Reichel and J. Esteve. E.W. acknowl- 
edges support from CNRS and Pohsh GRF: N N202 128539. 



A Some useful relations 

Here we collect useful commutation relations and recall the expression of the Bogoli- 
ubov quasi-particle annihilation operators c^-k after the mixing pulse in terms of the 
atomic annihilation and creation operators. 

The numerical coefficient in this formula is given for a cubic quantization box with 
periodic boundary conditions, see discussion in section 7.8 of [36| . 
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Fig. 8. (Color online). Condensate squeezing (blue upper curves) and squeezing of the 
total field (red middle curves) as a function of time. Full lines: classical field simulation. 
Dashed lines: Complete classical field Bogoliubov approximation implemented numerically. 
The two-mode result (black lowest dotted lines) is shown for comparison. Upper graph: 
kBT/pg = 7.83, {Nnc)/N = 0.05, N = 983040, pg = 693.Gh^ /mV^^^ , ^fp^ = 4.17 x 10"*. 
Lower graph: kaT I pg = 0.5, {Nnc)/N = 0.02, N = 2733750, pg = 13715^"^ /mV^^^ , 
= 1.32 X 10"^ 



Commutation relations: In our lattice model, the field operators obey discrete 
bosonic commutation relations: 



dV 



Vcr, a' ~ a,b 



(76) 



The hermitian condensate number and phase operators obey, for a, a' = a, b: 

[N^O,0a']=lSa,a' (77) 
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The fields of the non-condensed modes, orthogonal to the condensate wavefunction, 
obey 



[^/j„±{r),i'^,^{r )] = — - 



dV 



V 



(78) 



The non-condensed fields do not commute with the the total atom number but they 
commute with all the condensate operators: 



(79) 



The number-conserving operators Aa- obey the same commutation relations (|78l) as 
the non-condensed fields, e. g. 



[^(r),it(r')]^|^-l ya = a,b 



(80) 



but, contrarily to them, they commute with the total atom number operators in each 
component: 

Vcr, cr' = a, 6 (81) 



Their commutation relations with the condensate operators are, for a, a' = a, b: 



= 0, 



(82) 



Finally, from the relation e~*''''/(7Vao) e*^° = f{Nao - 1) resulting from ([771), for a 
generic function /, we have in the large N, and thus large Nao limit: 



1 



O 



f 1 



(83) 



Bogoliubov transformations: By projecting ((29)) over the plane waves, we obtain 
after the pulse (t > 0): 



(84) 
(85) 



where dk and 6k annihilate an atom with internal state a, b and wave vector k. The 
inverse relations are useful: 



Cak = e "flk (yfc - a'_y^e " Vk 



(86) 
(87) 



with Uk, Vk given by ([50)1 . Note that it is often convenient to express Uk and Vk in 
terms of Sk, so that for example 



1 



1 



U'k + Vl = -{si + 2) and 2UkVk = -(s^ - s,') 



(88) 
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B After-pulse values of the condensate phases and quasi-particle 
annihilation operators 

To determine the condensate phase operator for the internal state a at time i = 0+ 
in terms of operators at t = 0~, we use the fact resulting from (|9I10[) that ao(0"'") = 



a, 



(0) .(0) 







6[)°Vv^ and 6o(0^ 



^0 



&[, Vv^- Then the definition ^ leads to 



(89) 



where the upper sign (— ) is for ct = a and the lower sign (+) for a = b, the modulus 
operator of an operator X is 



\X\ = X^'X 



1/2 



and we have introduced the number conserving operator 6q ' = e 
panding in the large Nj^^ limit, we obtain : 



^a(0+)-^i°^ =Herm 



TV N^:> 



o 



f 1 



(90) 

By ex- 



(91) 



We have introduced the decomposition 6g ~ x + iy, where x and y are hcrmitian 
operators, the usual notation { , } for the anticommutator and the notation Hcrm X = 
{X + X'l^)/2 for the hermitian part of an operator X. This leads to ([57)) . 

To determine the quasi-particle annihilation operators Co-k at time t = 0+ in 
terms of operators at time i = 0~, we use the relations (|86I87|) to express them in 
terms of the atomic creation and annihilation operators at time t = 0^, that are in 
turn expressed in terms of their values at 0^ thanks to (|9ll0p . One also needs the 

expansion exp[i0<,(O+)] = exp[i^i''^] [1 + i(^^(0+) - ^i°^) -hO(l/A^)] deduced from ^ 



and the Hausdorf formula. With the short-hand notations a^^ 

m 



exp(— 



(0)n;-(0) 



and 



:exp(-*^i°))Sf: 



cak(0+) 



V2 



V2 



V2 



V2 



TV 



(92) 



where the upper, — sign is for a = a and the lower, -I- sign if for a — h. One then 
introduces the operators 



ik = and Bk = ^.iL"^-ViiLT 



(93) 



This directly gives (|40)) . Expressing a|^°'' and a^^^ in terms of the pre-pulse quasi- 
particle annihilation and creation operators thanks to the t = Q~ equivalent of ([S3]), 
gives (|39p . Restricting the accuracy of (|92p to 0(1) included, gives 



C Correlations of and 

Some useful properties of the operators defined in p9l40p [or cquivalently in ([M|] are 
given here. The commutation relations of these operators are bosonic. This means 
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that the only non-zero commutators (considered among all possible values of k) are 



= 1 



(94) 



To calculate averages of products of A and B operators (here at equal times) one 
can use the Wick theorem as the initial density operator for {c|j^^} and {6^°''} ^ 
Gaussian. All the non-zero correlations involving the can be deduced from 



1 



,(0) 



,(0) 



r„(o)2 


4 


2 ~ 


J0)2 


" J0)2 


4 




J0)2 



(95) 
(96) 



with s[°'' and Sk defined by (|28l30p . All the non-zero correlations involving the i?k 
can be deduced from 



'k^-2 



(i?kB-k) = -;7fci4 = --^4-^ 



(97) 
(98) 



All the crossed second moments, for example of the form (AB) or (A^B), are zero. 



D Double expansion of some operators 

As explained in the main text, to have a vanishing error on the squeezing parameter 
in the thermodynamic limit, it suffices to determine the operators Sy and Sz up 
to sa N^/^ included. 

Case of Sz- The operator Sz = {Na — Ni,)/2 is a constant of motion, it can be 
evaluated at t = 0+, and related with (|9ll0p to the fields at t = 0~ . Then one uses 
the modulus-phase representation for the condensate operator in a and one introduces 
the number-conserving fields and e''^-^ for the non-condensed modes, whose 
Fourier components can be expressed in terms of the operators A^ and Sk through 
The only approximation is then to neglect the commutator of y with (A'^o')^^^, 
which is 0{1/N^/^), to obtain 



Na^m-^^ U va1°)| - J2 [iui + )(4^k + A^Bi) 

+2UkVk{AlBi^ + A^B^k) (99) 



We recall than x and y are defined below (PT|) . The operator Sz has a zero expectation 
value, and a variance exactly equal to N/A, as already found by a more direct method 
in (1311), so we reach the estimate 



Sz « Ari/2 



(100) 
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With more lengthy calculations, we now deduce Sy from the antihermitian part 
of S-^- written in the form P5)) and we then obtain Sy and {Sy,Sz}- 

The phase difference: We first evaluate the scaling of the phase difference 9a — Ob 
in the thermodynamic limit from the writing psp . The contribution of the phase 
difference at time t = 0^ scales as according to (|37)) . The contribution pro- 

portional to Na — Nb scales in the same way, since the total number difference « N^^"^ 
for the binomial distribution after the 7r/2 pulse. The same conclusion holds for the 
contribution proportional to V, see (|T7|) . We reach the important conclusion that, for 
a finite time t in the thermodynamic limit, 

^^-^"^J^ (101) 

As Y + ^ is 0{N), it suffices to expand the exponential in (031) to first order included 
in the phase difference to obtain 



^4 



1 - ii§a -eb) + o 



N 

-+FR + iFj] (102) 



where we have split F ~ Fji + iFi in terms of the hermitian operators Fft. and Fj. 
The antihermitian part of F: The operator F/ directly contributes to the anti- 
hermitian part of 5+, so it has to be evaluated up to « N^/"^ included. Its exact 
expression is 

P^ = k^dV (it _ AlAa) ^^Y. (4c.k - h.c.) (103) 

This corresponds to a complex scalar product between the bicomponent fields {Aa^ A\) 
and (yife, Al). The Bogoliubov equations of motion for spin state a conserve this scalar 
product [3T] . Due to the a ~b symmetry, the coefficient of Bogoliubov equations of 
motion are the same for the two internal states, and Fj is a constant of motion 
within Bogoliubov theory. We can thus evaluate it at time 0"*", taking into account 
the corrections to c\^ci{^^) due to the small condensate phase change induced by the 
pulse, as in ([5^ : 

- k - BlA.) -\h ^^-^^ (104) 

where the operator y is defined below (j9ip . This correction involving y is important to 
ensure that (5^(0+)) — N/A as it should be. The operator Fj has a zero expectation 

value, this is why the same phenomenon as for the operator 2? occurs. Calculating its 
variance, which is dominated by the contribution of the sum over k in (|103|) . 



(105) 



we get as in (|47| the estimate 

Fi « (iVeBog)'/' (106) 
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The hermitian part of F: Contrarily to Fj, the operator alone cannot contribute 
to the antihermitian part of 5+, it has to be multiplied at least once by the phase 
difference operator. To obtain Sy up to « A^^/^ included, we thus need up to « iV 
included. To this end we decompose iV^o after the pulse, for <t = a or b, as follows: 

N,a = ^+5N„, 5N„={n,--\- N,^ (107) 



and we expand the square root in F (|44p in the large N limit to obtain 

: — iV 1 

^{Nao + l)Nbo = y + ^(1 + 5Na + SNi,) 



2 2 

{l + 5Na-5Nb? 



AN 



(108) 



In the second contribution in the right-hand side of pOSp , we can replace Na + 
with N . Since SNa — SNi, scales as iV^/^, the third contribution scales as N° and is 
thus negligible at the required order. The terms in the . . . are of too high order in 
Esizc or in eBog to be relevant. We conclude that, for our purposes, we can take 

F« ^ _i ^dF(it _ Al){Aa A,) (109) 



By expanding the fields Aa and Ab over the Bogoliubov modes, we obtain Fn in terms 
of the quasi- particle annihilation operators Co-k at time t > 0. Using we can relate 
these operators to their value at time 0"*", that we can replace by the leading order 
expression ([55)) to obtain 



Fr^-Y, [Vif + {Ul + 14')^k^k + UkVk (e-2-'^*/''^kS_k + h.c.)] (110) 



Its expectation value is 



{FR) = -Y,AUMsir,^u:kt (111) 

k#0 



with ojfe = efc/fi. Even if the _B's correspond to vacuum fluctuations, we still find 
(replacing the sum by an integral over K'^, which is convergent, and making the 
change of variable k = K/^heai, where ^hoai is the healing length) that {Fr)/N scales 
as (pa^)^/^, which is O(eBog)- For simplicity, we shall forget about this detail and 
consider that {Fr) w iVeeog- Using Wick's theorem we have also determined the 
variance of Fr, 

Var Fr^Y^ ^UlV^ sin^ Ukt (l + WlV^ sin^ LOkt) (112) 

Since (C/feVfc)'* diverges as for k 0, with a similar reasoning as for the discussion 
around ([751) . we find that Yav Fr is O {(NeBog)^^^) uniformly in time. We summarize 
these estimates by the writing 

Fr « NeBos ± (^eBog)'/^ (113) 



28 



Will be inserted by the editor 



As we said, we need to estimate Fn up to w iV included. This means that the fluctu- 
ations of Fn, that are N^^^ times smaller, are negligible and the operator Fn can be 
replaced by its mean value {Ff(} . 

Operators Sy, etc: From the antihermitian part of p02[) . and from the estimates 
(|101I1131 1106p . we can approximate Sy up to the terms w N^^^ included as 



Sy ^ Fj 



N 



(114) 



Squaring this expression, and neglecting terms of order larger than one in eeog, we 
finally obtain the expectation value {Sy)/N with an accuracy up to « e*^^^^ and « 
EBog included, see Taking the anticommutator of pi4p with Sz and then the 

expectation value gives (j53p . with the additional simplification that (FjSz) is purely 
imaginary and cancels out in the anticommutator [this results from (j99p and p04p . 
and in particular from (yx) = l/(4i)]. 

To conclude this Appendix, we give an expectation value useful for subsection [ 



E4 

k#0 



(115) 



where the A' A expectation value is given by We also give the expression of the 
difference of non-condensed atom number operators at t > 0: 



N, 



+ 2UkVk (ikS^ke 



-2iuJkt 



+ h.c. 



(116) 



and the corresponding variance written as a sum of non-negative terms, useful for the 
discussion below Ea.([7i|: 



YaT{Na^-Nbi_)^{N'^^l) + 




y- 

k^tO 

„(0)2\ 



,(0) 



- I sin Wfci 



,(0)2 



,(0)2 



COS Wfci 



(117) 



E With the oscillating terms in Qa — K 

As announced above ([55]) , we give here the analytical result for ^ (t) (with the double 
expansion technique) without performing the approximation used in the main text of 
the paper. The temporally oscillating terms in the phase difference operator are now 
kept, which amounts to replacing T) in (|35p with T)to\. = 2^ + ^^030 with the oscillating 
contribution 

23osc(t) = - E (e— '=*ikB_k + h.c.) (118) 

k#0 

where 0;^ = Cfc/fi- In the renormalized time (|59p one has also to replace T) with I?totj 
which involves the new expectation value 

({^.,Posc}) = E 4^^^ ((iki-k) + C/fcVl.) (119) 
k^o 2""=* ^ ' 
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that will thus be added to the contribution (jllSp [see (1^5]) for the A|j.^_k expectation 
value]. One can show that t{{S ztVosc}) is uniformly bounded in time, so it contributes 
to T as a time dependent small temporal shift. The result ((60|) is replaced by 



As expected, again, 2? was replaced by 2?tot- There is also an extra term, ({t), which 
is simply the value of the last term of ([5^ at time t minus its value at time 0+. This 
difference is no longer zero, because (j56|) no longer holds when V is replaced with 
T>tot, since I?tot has imaginary components. We find 

Ncit) = - Vi"' (4") + 1) (121) 

Note that C(i) is uniformly bounded in time, as (Ffl)/iV, it is thus not particularly 
relevant. 

More significant deviations may come from the occurrence of {'D'^^^) that differs 
from the original (P^) by the two terms 



k^^O 

j2 



LOkt 

+Ul - 2UuVu (^k^-k) cos 2L0kt (123) 



These two terms however are 0(eBog/'''^) in the long time limit. At the relevant times 

1/2 

r « l/^Bog! "where the multimode nature of the field starts limiting the squeezing, 
their contributions to ^^ot a-i'S 0[e^^^) and negligible. 

To summarize, the inclusion of the non-oscillating terms in the phase difference 
operator does not change at all the long time limit of (which, importantly, is its 
infimum in the Bogoliubov approximation): Eq. (j6ip is unchanged. At intermediate 
times, it gives small deviations. For the extreme case ksT/pg = 10 and {pa^y^^ = 
10~^, where the non-condensed fraction reaches 10%, we find a maximal relative 
deviation of 2% between ^^(t) of ((60|) and the more accurate ^tot(Oi ^■t a time pgt/h ~ 
1.5 when is still a factor ~ 4 above its minimal value ~ 0.1. 
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